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We consider a dilute gas of granular material inside a box, 
kept in a stationary state by shaking. A wall separates the 
box into two identical compartments, save for a small hole 
at some finite height h. As the gas is cooled, a second order 
phase transition occurs, in which the particles preferentially 
occupy one side of the box. We develop a quantitative theory 
of this clustering phenomenon and find good agreement with 
numerical simulations. 
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One of the most outstanding features of a gas of granu- 
lar material is its tendency to spontaneously form highly 
concentrated regions or clusters ||l|-|[. So even in its 
gaseous state it behaves fundamentally different from a 
molecular gas, which keeps its uniform density. Apart 
from throwing light on the noncquilibrium properties of 
a granular gas, understanding this clustering instability 
is of major technological importance. Imagine a flow of 
rocks down a chute: whenever a very dense region has 
formed due to the instability, the rocks may easily get 
entangled and the flow is stuck. 

The clustering instability is caused by the distinguish- 
ing feature of a granular gas in contrast to a molecular 
gas, namely that a fraction of the kinetic energy is ir- 
retrievably lost upon a collision, and simply heats the 
particles. In this paper, we will concentrate on the sim- 
plest case of a dilute gas for which collisions are binary, 
so the rate of energy loss grows quadratically with the 
density. This means a dense region rapidly cools, in- 
creasing the density even more according to the equation 
of state. Goldhirsch and Zanetti ||] used a hydrodynamic 
description of a dilute gas or "rapid granular flow" jj-g] 
to show that this mechanism leads to a long-wavelength 
instability of a homogeneous assembly of inelastic par- 
ticles. However, in the absence of driving, the collapse 
quickly becomes complete and the original assumptions 
break down. In the system described in this letter, a 
stationary state is maintained by external driving. It is 
thus simple enough to allow for an analytical descrip- 
tion, while it keeps the central aspects of the clustering 
phenomenon. 

The experiment, first described by Schlichting and 
Nordmeier M, consists of a box of base area 12cm^ 
and height 20cm, mounted on a shaker, and filled with 
N = 100 plastic particles of radius r — 1mm (see Fig. ^). 
The box is separated into two equal parts by a wall which 
has a narrow horizontal slit at a height h = 2.3cm. When 
the shaker is operating at full power, the amplitude is vi- 



bration is approximately A = 0.3cm and the frequency 
/ = 50Hz. Even if all particles are on one side initially, 
they immediately distribute equally to both sides. Low- 
ering the frequency below a critical value of 30Hz the 
symmetry is spontaneously broken, and particles settle 
preferentially on one side. Evidently, the wall greatly en- 
hances the clustering described above, since it prevents 
a direct exchange of particles which would break up the 
cluster. The grains thus act as Maxwell's demon pO[ , 
who preferentially lets particles pass from left to right or 
vice versa. As a result, a more ordered state is formed in 
which most particles are on one side. The demon must 
then absorb entropy, a role which in our system is as- 
sumed by the sand grains. Still another interpretation 
would be that of a dissipative structure, which in the 
stationary state is maintained by a flux of entropy pd| . 
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FIG. 1. A drawing of the experimental setup. The two 
sides of the box are connected by a hole at height h. The 
picture is taken below the bifurcation, so most particles have 
settled on the right hand side. As a result, the gas sinks to 
the bottom, reducing the flux. 

The condition of stationarity is that the total flux be- 
tween the two compartments is zero, i.e. the fluxes going 
from one side to the other cancel: 



Fi^rih) = Fr^eih). 



(1) 



In the limit that the connecting hole is small, it remains 
to find the density and temperature distribution in a 



one-dimensional column of a vibrated granular assembly, 
from which we calculate the flux of particles leaving a 
compartment at the height of the hole h. Below we show 
that (H) has a single solution for strong driving, but two 
asymmetric solutions become stable as the driving is low- 
ered. 

The problem of a shaken granular gas in a gravita- 
tional field has recently been treated in a number of pa- 
pers p^-p^ . Equations of motion for the number density 
n(z), pressure p(z), and granular temperature T(z) can 
be found from conventional kinetic theory EIJ^ . The posi- 
tion z is measured from the bottom of the container and 
the granular temperature is defined a,s T =< v^ > /d, 
where d is the dimension of space and v the velocity of 
a particle. This definition of temperature is customary 
for granular media, but differs from the usual molecu- 
lar temperature, which is recovered by formally setting 
ks = rn. Now the stationary equations become in the 
dilute limit 



VbT{0) + ,/2/^vlT{0y/^ 



(3) 



p — mnT, dzP = —mgn, 



(2) 



The first equation is the usual equation of state, the sec- 
ond the force balance, where —mg is the force on a single 
particle. The third equation is the balance of heat flux 
and dissipation due to inelastic collisions. As usual, the 
thermal conductivity kT^'^ is proportional to the aver- 
age particle velocity. In the simplest possible model of 
hard and smooth particles the energy loss in one collision 
is proportional to (1 — e'^)T on the average. The coeffi- 
cient of restitution e measures the reduction v'„ = —evn 
in the normal velocity of the particles. Together with the 
number of collisions being proportional to ri^T^/'^ this ex- 
plains the form of the loss term on the right hand side of 
the third equation (||). 

To allow for a better comparison with numerical sim- 
ulations, we will consider the two-dimensional case of 
circular discs of radius r, for which the coefficients are 
found to be g K = Tr^^/^m/r, D — 47r^/^mr(l — e). 
Next, we have to supply boundary conditions. To mini- 
mize wall effects, which are not essential to our problem, 
all wall collisions are assumed to be elastic. The top is 
left open. For simplicity, the bottom of the container is 
taken to move in a sawtooth manner, such that a collid- 
ing particle always finds it to move upward with velocity 
Vb = Af . Finally, the amplitude A of the vibration is 
assumed to be very small compared with the mean free 
path, so that the bottom is effectively stationary. In sum- 
mary this means that the z-component of the velocity of 
a particle colliding with the bottom is changed accord- 
ing to v'z — 2wfc — Vz- The approximations of (||) imply 
that the velocity is close to a Maxwell-Boltzmann distri- 
bution, which allows to calculate the rate of energy input 
per unit width (or unit area in 3D) to 



We will see below that the first term in (y) dominates, 
the second typically being smaller by a factor of 10 in 
our simulations. Thus the two boundary conditions for 
z — Q become 



p(0) = -gmN, Q = -KT^/^dzT{0). 



(4) 



The first equation comes from integrating the force bal- 
ance, and N is the total number of particles per unit 
width of the box. The second equation balances the en- 
ergy input with the heat flux out of the bottom according 
to®. 
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FIG. 2. Temperature and density profiles of A^ = 320 
grains in a box of width w = 1.6. The particles are circular 
disks of radius r = 0.01 and coeflicient of restitution e = 0.95. 
The velocity of the bottom is v-t = 0.149, acceleration of grav- 
ity g and mass m is normalized to one The full line is the result 
of a particle simulation, the dotted line is the present theory, 
and the dashed line is (0). 

Since (0) is of third order, another boundary condition 
is needed for a unique solution. The missing third con- 
dition is found by observing that the solution of (pi) for 
large z is of the form T'^/^(z) = cz + 6, so c must be zero 
to prevent T from becoming negative or infinity. This is 
the desired third condition, and solutions are easily found 



by shooting for a temperature profile which is asymptot- 
ically constant. The resulting temperature and density 
profiles are shown as dotted lines in Fig. g for a typical 
set of parameters and are compared with a numerical sim- 
ulation (solid lines). The simulation uses an event-driven 
code |15i for smooth and hard particles, with the normal 
velocity reduced by a factor of e as described above. The 
agreement is quite good, but gets worse if the number of 
particles is reduced, and deteriorates even more in three 
dimensions. The reason is that even for the parameters 
of Fig. the temperature changes significantly over the 
length of the mean free path. As a result, the distribution 
of the Wz-velocities of particles hitting the bottom devi- 
ates significantly from a Maxwell-Boltzmann distribution 

at r(o). 

The solution given in ||13| assumes that the tempera- 
ture is constant, which is asymptotically correct for very 
small inelasticity. The resulting profiles are the dashed 
lines in Fig. ^ Evidently, the present theory is more 
accurate, but the constant temperature solution has the 
great merit of simplicity and still contains the essentials. 
Namely, at constant temperature the density is an ex- 
ponential and Too is found from a balance of the energy 
input and dissipation: 
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(5) 



Here we have neglected the second term of (|^) in favor 
of the first, which is again consistent for e « 1. 

Returning to the original stability problem, the flux 
through a hole of area S is found to be F — 
Sn{h)^jT{h)/2Ti. Using this formula, and solutions of 
(^-(0) to find the profiles, one can look for solutions 
of (P^. Since the total number of particles N is con- 
served, these must be sought subject to the constraint 
N = Ni + Nr, where the overbar denotes the number of 
particles relative to the half-width of the container. 

Figure S shows the solutions of (0) as a function of 
h for N = 225 as the solid line. The instability fol- 
lows a typical pitchfork bifurcation, so the average of the 
asymmetry parameter e = {N^/r — N)/N increases con- 
tinuously when h is raised above a critical value. The 
circles, which are the results of a numerical simulation 
with N = 360 particles averaged over time, show quite 
satisfactory agreement. Figure Q illustrates the temporal 
behavior of the asymmetry parameter e. The symmetric 
state rapidly becomes unstable and the system fluctu- 
ates around its new equilibrium position. Since there 
are more collisions on the denser side the temperature is 
lowered, causing particles to sink to the bottom. Thus 
a non-symmetric state becomes possible because the flux 
F is no longer a monotonously increasing function of the 
number of particles as it would be in equilibrium. An- 
other hallmark of a nonequilibrium system is that the 
stationary state is described by a flux balance (|^) , while 



the temperature on either side of the hole is not equal, 
as it would be in thermal equilibrium. 
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FIG. 3. The bifurcation of the asymmetry e as function of 
the height h. The full line is the result of the present theory, 
circles are numerical simulations with A^ — 360 particles in a 
box of half-width 1.6, averaged over time. All other parame- 
ters are the same as in Fig. ti. The dotted line is the result 
(I?!) of the simplified theory. 
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FIG. 4. The temporal evolution of the asymmetry param- 
eter e corresponding to the data point h — 0.9 in FigJll, but 
with only N — 180 particles. 

Further insight can be gained from the simplified solu- 
tion (||). A straightforward calculation leads to 

dte = FqN 

(e-l/2)2e-^(^-i/2)' 



(e + l/2)2e-''(^+i/2)^l ^^^ (g) 



where ^ is a noise term to be considered later and Fq 
is a constant. The stationary solutions for small noise 
are determined by the vanishing of the angular brackets, 
which is controlled by a single parameter 



H = ATTghr^{l~eyNyvt. (7) 

If /i > 4, the equation [...] = has three roots, and just 
above the bifurcation the asymmetric solutions are de- 
scribed by < e >= ±y^3(^i^4)/16, which is included as 



the dashed hne in Fig. |^. While there is an offset between [4] J. T. Jenkins and S. B. Savage, J. Fluid Mech. 130, 187 

the simplified theory and the simulation, it does describe (1983). 

the form of the bifurcation fairly well using the single [5] J- T. Jenkins and M. W. Richman, Phys. Fluids 28, 3485, 

parameter /x. It is intuitively clear that by increasing the (1985). 

number of particles, the inelasticity, or the height of the I^] J. Jenkins and M. Richman, J. Fluid Mech. 171, 53 



wall separating the two compartments, phase separation 



(1986). 



, , ^? ^. .^, , . ' f^ . , J ,. [7] C. S. Campbell, Annu. Rev. Fluid Mech. 22, 57 (1990). 

IS enhanced. Heatmg, as pomtcd out m the mtroduction, jgj ^ g^^^ ^^^ ^ Goldhirsch, J. Fluid Mech. 361, 41 (1998). 

favors a symmetric state. jgj ^ j gchlichting and V. Nordmeier, Mathematisch 

Finally, it is possible to include fluctuations in our de- Naturwissenschaftlicher Unterricht 49, 323 (1996). 

scription, which are important if the number of particles [lo] A. S. Leff and A. F. Rex, Maxwell's demon: entropy, 

is just a few hundred as in Fig. ||. In particular this information, computing, Adam Hilger, Bristol (1990). 

leads to a softening of the transition, since immediately [11] P. Glansdorff and I. Prigogine Thermodynamic theory 

below the transition the system can switch between the of structure, stability and fluctuations, Wiley, London 

two possible states. Assuming that the particles passing (1971). 

through the hole are uncorrelated, the noise ^ in (§) is [12] J- Lee, Physica D 219, 305 (1995). 

Gaussian white noise pl on a coarse-grained time scale, I^^] V. Kumaran, Phys. Rev. E 57, 5660 (1998). 

, il_li t> j^^j g McNamara and S. Luding, Phys. Rev. E 58, 813 

^"^ (1998). 

<- ^(t)^(t') >= Fo [15] D. C. Rapaport, J. Comp. Phys. 34, 184 (1980). 

, , ,„,2 „ , , ,„,2n fl6] H. Risken, The Fokker-Planck equation. Springer, Berlin 

(e-l/2)V''(-V2)^ + (e + l/2)V^(^+i/^)^]<5(t-t'). (8)^ \i984). 

Note that the constant Fq can be eliminated by rescal- 
ing time, so the strength of the noise is controlled by the 
total number of particles alone. By considering the fluc- 
tuations around the local minimum to Gaussian order, we 
find ^ ^< (e- <e>)2 > = [AN{p. - A)]-^/"^ . Adjust- 
ing this formula to the critical value of fXcr = 5.4 found 
from simulation this gives yj< (e— < e >)^ > = 0.045 
in reasonable agreement with Fig. ^. Of course, equa- 
tions (^-(H) also allow for a more detailed analysis of 
transitions between the two asymmetric states and many 
more questions relating to the critical fluctuations near 
the transition. This will be considered in more detail in 
future publications. 

In conclusion, hydrodynamic equations are a very valu- 
able tool to describe a dilute granular gas. The greatest 
problem lies in the formulation of the boundary condi- 
tions. The clustering instability can be understood in 
terms of a very simple experiment, which we model in a 
static, one-dimensional description. 
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